In [1]:
%pylab inline
import matplotlib.pyplot as plt
Populating the interactive namespace from numpy and matplotlib

Polariton, phonon-photon interaction

C. Kittel: Quantum Theory of Solids, Chapter 3, page 43, Eq. 39.

Dispersion relation:

$\omega = \sqrt{\frac{1}{2\varepsilon_{\infty}}(\omega_t^2 \varepsilon_0+ c^2 k^2)\pm \Omega^2}$, where $\Omega^2 =\sqrt{\frac{1}{4 \varepsilon_\infty^2}{(\omega_t^2 \varepsilon_0+ c^2 k^2)}^2 - \omega_t^2\frac{c^2 k^2}{\varepsilon_\infty}}$

Paraméterek:

e0 ---> $\varepsilon_0$ dielektromos állandó $\omega = 0$ frekvencián

ev ---> $\varepsilon_{\infty}$ dielektromos állandó $\omega = \infty$ frekvencián

omt ---> $\omega_t$ tranzverzális otikai módus frekvenciája

oml ---> $\omega_l$ longitudinális optikai módus frekvenciája

In [2]:
# Az abra kimentesehez az alabbiakat a plt.show()  ele kell tenni!!! 

#savefig('fig_rainbow_p2_1ray.pdf');  # Abra kimentese
#savefig('fig_rainbow_p2_1ray.eps');  # Abra kimentese

# Abra es fontmeretek
xfig_meret= 9   #    12 nagy abrahoz
yfig_meret= 6    #   12 nagy abrahoz
xyticks_meret= 15  #  20 nagy abrahoz
xylabel_meret= 21  #  30 nagy abrahoz
legend_meret= 21   #  30 nagy abrahoz
In [3]:
### C. Kittel: Quantum Theory of Solids, Chapter 3, page 43, Eq. 39.  

def omega(k,params):
    
    e0,ev,omt = params
    
    tmp=sqrt((omt**2*e0+k**2)**2/4/ev**2-omt**2*k**2/ev)
    om1=sqrt(1/2/ev*(omt**2*e0+k**2)-tmp)
    om2=sqrt(1/2/ev*(omt**2*e0+k**2)+tmp)
    return ([om1,om2])

### epsilon(omega) dielectric constant with damping
### Ashcroft-Mermin: Solid Sate Physics, page 547, Eq. (27.57)

def epsilon_om(om,params,gammas):
    e0,ev,omt = params
    gamma, = gammas

    tmp = ev + (ev-e0)/((om/omt)**2-1+1j*(om/omt)*gamma)
    return tmp
In [4]:
e0=4;
ev=2;
#oml=4;
omt=3;
omlcalc=sqrt(e0/ev)*omt;
print("omega_l = ", omlcalc)

params=[e0,ev,omt]

print(params)
[o1,o2] = omega(10.,params)
o1,o2
omega_l =  4.24264068712
[4, 2, 3]
Out[4]:
(2.7256814723711122, 7.7827154972524113)
In [5]:
gamma =0.2
gammas=[gamma]
eps = epsilon_om(0.7,params,gammas)
eps
Out[5]:
(4.1100190780547265+0.10413725179588545j)
In [6]:
k = linspace(0,7,10)
o1, o2 = omega(k,params)
o2
Out[6]:
array([ 4.24264069,  4.26057568,  4.31569762,  4.41174033,  4.55404953,
        4.7479234 ,  4.99630797,  5.29799685,  5.64763602,  6.03744483])

A C. Kittel könyvben a 43. oldalon a (41) képlet rossz, helyesen: $\omega = ck/\sqrt{\varepsilon_{\infty}}$

In [7]:
Npoints=100;
kmax=17;
k = linspace(0,kmax,Npoints)

omphoton=k/sqrt(ev)   

omlplot=linspace(omlcalc,omlcalc,Npoints)
omtplot = linspace(omt,omt,Npoints)
o1, o2 = omega(k,params)

figsize(xfig_meret,yfig_meret)

plot(k,o1,color='b',label='$\omega_{1}(k)$')
plot(k,o2,color='r',label='$\omega_{2}(k)$')
plot(k,omphoton,ls='--',color='k')
plot(k,k,ls='--',color='g')

plot(k,omlplot,ls='--',color='k')
plot(k,omtplot,ls='--',color='k')

legend(loc='upper left',fontsize=legend_meret)


annotate(r'$\omega_{l}$', xy=(.7,2.), xytext=(-1.2,omlcalc),fontsize=21)
annotate(r'$\omega_{t}$', xy=(.7,2.), xytext=(-1.2,omt),fontsize=21)

xlabel(r'$c k$',fontsize=20)
ylabel(r'$\omega_{1,2}$',fontsize=20)
xlim(0,kmax)
ylim(0,10);
#xticks(fontsize=10)
#yticks(fontsize=15);


grid();
In [8]:
### Ashcroft-Mermin: Solid Sate Physics, page 551, Fig. 27.6
###  epsilon(omega) dielectric constant with NO damping

Npoints=100;
ommax=17;
x = linspace(0,ommax,Npoints)

epsv= linspace(ev,ev,Npoints)


gamma = 0.0
gammas=[gamma]

subplot(1,1,1)

epsreal = epsilon_om(x,params,gammas).real

plot(x,epsreal,label=r'$\mathrm{Re}(\epsilon(\omega))$')
plot(x,epsv,ls='--',color='k')

legend(loc='upper right',fontsize=legend_meret)
cim=r'$\epsilon(\omega), \mathrm{no\,\, damping}$'
title(cim,fontsize=legend_meret)

annotate(r'$\epsilon_{0}$', xy=(.7,2.), xytext=(-1.2,e0),fontsize=21)
annotate(r'$\epsilon_{\infty}$', xy=(.7,2.), xytext=(-1.2,ev),fontsize=21)


xlim(0,15)
ylim(-2,10)
grid();

print("gamma = ",gamma)
gamma =  0.0
In [9]:
### Ashcroft-Mermin: Solid Sate Physics, page 551, Fig. 27.6
###  epsilon(omega) dielectric constant with damping

Npoints=100;
ommax=17;
x = linspace(0,ommax,Npoints)

epsv= linspace(ev,ev,Npoints)


gamma =0.1
gammas=[gamma]

cim='Damping'

subplot(1,2,1)

epsreal = epsilon_om(x,params,gammas).real

plot(x,epsreal,label=r'$\mathrm{Re}[\epsilon(\omega)]$')
plot(x,epsv,ls='--',color='k')

annotate(r'$\epsilon_{0}$', xy=(.7,2.), xytext=(-2.7,e0),fontsize=21)
annotate(r'$\epsilon_{\infty}$', xy=(.7,2.), xytext=(-2.7,ev),fontsize=21)

legend(loc='upper right',fontsize=legend_meret)

title(cim,fontsize=15)

xlabel(r'$\omega$',fontsize=20)

subplot(1,2,2)

epsimag = epsilon_om(x,params,gammas).imag

plot(x,epsimag,label=r'$\mathrm{Im}[\epsilon(\omega)]$')
plot(x,epsv,ls='--',color='k')

legend(loc='upper right',fontsize=legend_meret)

title(cim,fontsize=15)

xlabel(r'$\omega$',fontsize=20)


xlim(0,15)
ylim(-2,10)
grid();
print("gamma = ",gamma)
gamma =  0.1
In [10]:
### Ashcroft-Mermin: Solid Sate Physics, page 551, Fig. 27.6
## Dispersion relation with NO damping

Npoints=100;
ommax=17;

gamma = 0.0
gammas=[gamma]

x = linspace(0,ommax,Npoints)

Np=20
kp=linspace(0,max(k),Np)
omlplot=linspace(omlcalc,omlcalc,Np)
omtplot=linspace(omt,omt,Np)

subplot(1,1,1)

kr = (x*sqrt(epsilon_om(x,params,gammas))).real

plot(kr,x,label='$\omega(k)$')
plot(kp,omlplot,ls='--',color='k')
plot(kp,omtplot,ls='--',color='k')

legend(loc='upper right',fontsize=legend_meret)
title("No damping, Real part of k")

annotate(r'$\omega_{l}$', xy=(.7,2.), xytext=(-1.7,omlcalc),fontsize=21)
annotate(r'$\omega_{t}$', xy=(.7,2.), xytext=(-1.7,omt),fontsize=21)

xlim(0,16)
ylim(0,11)

grid();
print("gamma = ",gamma)
gamma =  0.0
In [11]:
### Ashcroft-Mermin: Solid Sate Physics, page 551, Fig. 27.6
## Dispersion relation with damping

Npoints=100;
ommax=17;

gamma = 0.1
gammas=[gamma]

x = linspace(0,ommax,Npoints)

Np=20
kp=linspace(0,max(k),Np)
omlplot=linspace(omlcalc,omlcalc,Np)
omtplot=linspace(omt,omt,Np)

subplot(1,2,1)

kr = (x*sqrt(epsilon_om(x,params,gammas))).real

plot(kr,x,label='$\omega(k)$')
plot(kp,omlplot,ls='--',color='k')
plot(kp,omtplot,ls='--',color='k')

legend(loc='upper right',fontsize=legend_meret)
title("Real part of k")

xlabel(r'$\mathrm{Re}[c k]$',fontsize=20)

annotate(r'$\omega_{l}$', xy=(.7,2.), xytext=(-2.7,omlcalc),fontsize=21)
annotate(r'$\omega_{t}$', xy=(.7,2.), xytext=(-2.7,omt),fontsize=21)

xlim(0,16)
ylim(0,7)

subplot(1,2,2)

#ki = (x*sqrt(epsilon_om(x,params,gammas))).imag
ki1 = (epsilon_om(x,params,gammas)).imag
ki2=x*sqrt(ki1)

plot(ki2,x,label='$\omega(k)$')
plot(kp,omlplot,ls='--',color='k')
plot(kp,omtplot,ls='--',color='k')

legend(loc='upper right',fontsize=legend_meret)
title("Imaginary part of k")

xlabel(r'$\mathrm{Im}[c k]$',fontsize=20)

annotate(r'$\omega_{l}$', xy=(.7,2.), xytext=(-2.7,omlcalc),fontsize=21)
annotate(r'$\omega_{t}$', xy=(.7,2.), xytext=(-2.7,omt),fontsize=21)


xlim(0,16)
ylim(0,7)

grid();
print("gamma = ",gamma)
gamma =  0.1
In [ ]: